New ab initio approach for high pressure systems with application to a 
new high-pressure phase for boron: perturbative momentum-space 

potentials. 
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Through the use of perturbation theory, in this work we develop a method which aUows for a 
substantial reduction in the size of the plane-wave basis used in density-functional calculations. 
This method may be used for both pseudopotentials and all-electron calculations and is par- 
ticularly beneficial in the latter case. In all cases, the approach has the advantage of allowing 
accurate predictions of transferability errors for any environment. Finally, this method can be 
easily implemented into conjugate gradient techniques and it is therefore computationally effi- 
cient. In this work, we apply this method to study high pressure phases of boron. We find that 
boron undergoes a phase transition from the ai2-B structure to the Qga-B structure, both of 
which are semiconducting. The a^a-B structure has lower energy than traditional mono-atomic 
structures, which supports the assertion that the metallic, and hence superconducting phase, 
for boron is much more complicated than a simple mono-atomic crystal. 
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I INTRODUCTION 

Experimental techniques are now able to probe con- 
densed matter systems at higher and higher pressures 
through the usCp.of diamond-anvil cells or dynamical 
shock methods.ElEI At these newly attainable pressures, 
structural and electronic phase transitions can occur, 
opening the door to exploration for new physical phe- 
nomena. 

Such systems offer an exciting avenue for first princi- 
ple calculations, as their predictions can not only follow 
but also sometinL^Jbaii-Desults from new experimental 
techniques .□'□Ef&Q'ErB&Ey When applying traditional 
first principle calculations to such systems, care must be 
taken, as most basis sets take advantage of the distinc- 
tion between the core regions and valence regions which is 
prominent at ambient conditions. However, when study- 
ing systems over a wide range of pressures, this distinc- 
tion vanishes, and both regions need to be treated on 
an equal footing. The following section discusses some 
problems that can occur when applying some of the most 
popular basis sets to high pressure systems. 

To overcome the above potential pitfalls and to al- 
low for accurate description of both the core and va- 
lence regions while keeping the mathematical stability 
and systematic convergence of the plane-wave basis, we 
have developed a method based on perturbation theory. 
This method can be used for direct all-electron calcu- 
lations, dramatically reduces the size of the plane-wave 
basis while properly accounting for all higher momen- 
tum states, accounts accurately for the core electrons and 



their interaction with the valence region, and automati- 
cally generates pseudopotentials which change with the 
environment thereby enhancing transferability. More- 
over, the effect of pseudizing in the crystal can be quan- 
tified; hence, accurate predictions for transferability can 
be made in any environment. Finally, the method can be 
implemented into current state of the art minimization 
techniques, such as the conjugate gradient method, and 
is therefore computationally efficient. 

As an application, we study high pressure phases of 
boron. Boron crystallizes into many complex structpjes- 
which are governed by the regular icosahedron.EilE3 
Boron is semiconducting at ambient conditions, turns 
metallic under pressure (« 170 GPa) and superconduct- 
ing-above 160 GPa.0 Theoretically, Mailhiot and cowork- 
eraJ were the first to apply density-functional techniques 
to the study of the phases oLhoron at high pressures. 
To date, however, such studiesO't3 have been limited to a 
modest selection of phases and were based on traditional 
techniques which artificially separate the physics in the 
core and valence regions. In this work, we apply our new, 
unbiased technique to a wider range of phases and find 
a new phase, the a^o-B structure, to be lower in energy 
at high pressure than any phase reported previously and 
to have important potential implications for the observed 
semiconductor-metallic phase transition associated with 
the superconductivity. 

We proceed as follows: Section |l| describes the various 
methods that can be used to study systems under high 
pressure, specifically, the pros and cons of each method. 
Then in Sections III and IV our method is developed. 
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Computational details are presented in Section [V|, Fi- 
nally, we apply our approach to study the high pressure 
phases of boron in Section W^. 



II METHODS 

Within current state of the art electronic structure cal- 
culations, various approaches can be taken, mainly differ- 
entiated by the choice in basis set. Each set has its own 
benefits and disadvantages. Seme of the rrpst common 
basis sets-ajce. the plane- waves ,t3 Gaussian,t3 linearized 
methods,BE3 and waveletsO We will now go over the 
various basis sets and focus on problems that may occur 
when each is applied to high pressure systems. 

The Gaussian basis set expands the wave functions 
in terms of linear combinations of Gaussians designed 
to represent the atomic orbitals well. Such a basis set 
has the advantageous property of expanding the wave 
functions in a small number of computationally efficient 
basis functions. There are, however, some drawbacks. 
Notably, such a basis set cannot be improved systemat- 
ically. This does not tend to be a practical problem for 
normal systems; however, such a problem can manifest 
itself for a system studied at various volumes, as the ba- 
sis set will span a different percentage of real space for 
systems whose volume differ appreciably. Because it is 
unclear of how to convergence such a basis set systemati- 
cally, it is unknown whether standard Gaussian basis sets 
will perform well at high-pressures, particularly because 
such bases are biased toward orbitals constructed at am- 
bient conditions. It is thus unclear if such a bias can 
hamper calculations at high pressures, where unknown 
phenomena may occur. 

Linearized methods generally fall under one of twa 
methods: the Linear Muffin Tin Orbital (LMTO)E£l 
method j-Qr the Linearized Augmented Plane- Wave 
(LAPW)El method. Of these two, the LAPW is more 
accurate and we will therefore concentrate on this ap- 
proach. When studying systems under high pressure, the 
most notable problem that can occur with the LAPW 
method is the treatment of the core electrons. At high 
pressures, certain core electrons need to be promoted and 
treated as valence electrons, and care must be taken with 
this method when promoting such electrons. Two widely 
used approaches to propapte such electrons are the mul- 
tiple wiBjdow approachlij and the localized orbital ap- 
proach.ll3 Both of these approaches could be particularly 
problematic for high pressure systems, particularly for 
first row elements. The multiple window method is not 
guaranteed to give a good solution,!!^ as different basis 
sets are used to expand the core and valence electrons. 
In the localized orbital approach, this problem is resolved 
by using one energy window. However, problems may ex- 
ist, as there are very large energy separations between the 
valence bands and the core bands, particularly for first 
row elements. 



The plane-wave approach, when used with the pure 
Coulomb potential, has the advantageous property that 
convergence of the energy can be systematically improved 
in a stable and controllable framework. Moreover, such 
a basis is unbiased and can therefore represent core and 
valence regions on an equal footing. However, such an 
approach has onliz_.hfiaiL. used for thej_aimplest systems, 
such as hydrogenllfl'LJ'EHl and lithiumEy as it requires a 
huge number of plane-waves to properly account for the 
singular Coulomb potential and the exclusion principle. 
Therefore such an approach is not practical beyond the 
simplest systems. 

The pseudopotential approach had been developed 
to cirpuniijeiLt the above difficulties with plane waves 
bases.Eilcj'Ej However, pseudopotentials are poorly 
suited for studying high pressure systems, iiecause, al- 
though some key properties are known ,Ej'Ej'E£I the trans- 
ferability of pseudopotentials is not systematically under- 
stood, particularly when the pscudized regions occupy a 
large percentage of the available volume as they do at 
high densities. More fundamentally, apart from nonlin- 
ear |Core corrections in the exchange-correlation poten- 
tial,E3 pseudopotentials completely ignore the core re- 
gion, which have i^ff n shown to play a vital role for high 
pressure systems Bq 

One of the most successful techniques, which keeps the 
benefits of the pseudopotential approach and can allow 
for core electrons and can recover details of the true wave 
function in tha.core region, is the Projected Augmented 
Wave (PAW)t3 method. This- method has been used 
to study Lithium and Sodiunij'Q under high pressure. 
Such a method is quite elegant; however, it must be con- 
structed from a referenced atomic state. Moreover, the 
construction of such a potential still requires a real space 
cutoff. In order to effectively study high pressure sys- 
tems, such a cutoff should be made quite small, so that 
the pscudized regions do not overlap. The approach then 
looses many of its computational benefits. 

The final type of basis set of which we are aware ia. 
the study of solid state systems are wavelet bases.EZI 
Such bases have the attractive property of using multi- 
resolution analysis to provide the necessary resolution in 
each region of space in a systematic and mathematically 
stable manner and thereby efficiently handle valence and 
core electrons exactly. However, such a basis set has just 
come to the fore and has not yet been adopted by many 
groups. 

We now describe a technique which has the advan- 
tage of using a plane-wave basis without pseudopoten- 
tials. This method has the added benefit of greatly re- 
ducing the size of the plane- wave basis set by generating 
the important high momentum states through perturba- 
tion theory. Below, we show that a properly constructed 
perturbation expansion can be quite accurate in this con- 
text. Finally, we show that this method can easily be im- 
plemented in current conjugate gradient techniques and 
therefore can be used with highly optimized minimization 
techniques. 
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Ill PERTURBATIVE DERIVATION 

In the next two sections the method will be derived. 
To gain insight into this method, it will first be derived 
for the simple case of directly diagonalizing the single- 
particle Hamiltoniaii. Then, connections will be made to 
previous approachEa In Section IV, the procedure will be 



implemented into the conjugate gradient framework, so 
that optimized calculations can be performed. Technical 
details will be discussed in Section For simplicity in 
the formal developments below, we sample the Brillouin 
zone at the Gamma point. Adding in k-point dependence 
is straightforward. 



A Derivation for fixed potentials 



Standard electronic structure calculations seek for the 
minimum of an energy functional £'[{C}] as a function 
of a set of basis function coefficients {C} subject to an 
orthonormality constraint. The variational derivative of 
this energy functional and constraint then lead to the 
standard eigenvalue problem 



HC = eC, 



(1) 



where H is the single-particle Hamiltonian (fixed and 
not updated self-consistently for the discussion in this 
section), C is an eigenvector and e is the corresponding 
eigenvalue. Note that we here assume that the basis set 
is orthonormal as our objective here is to work with plane 
waves. 

Separating the coefficients into two groups, those de- 
scribing behavior in the low spatial frequency space P 
as and those describing behavior in the high spatial 
frequency space Q as C"^, also separates the couplings 
which H describes into four groups, couplings from low- 
frequency to low- frequency H^^, from low- to high- fre- 
quency H^^ ^ from high to low H^^ and from high to 
high iJ'3Q. The above equation then decomposes simply 
into 



(2) 
(3) 



Solving for the Q-space components of the eigenfunction 
in terms of the P-space components leads to 
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(4) 



where the fraction mcans-tbej inverse of an operator. 

A standard techniqucEj'E3'L3 is to substitute Equa- 
tion (U) into Equation in order to generate an effective 
eigenvalue equation for the P-space, 
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HQQ 



HQP)cP = eC' 



(5) 



This equation reduces the problem to the P-space only. 
If standard diagonalization techniques are used, this re- 
duces the time to diagonalize the Hamiltonian by a factor 
of ( jv-P^jvg )'^; where iV" is the number of basis functions 
in the a-space. However, it is well known that such a de- 
composition does not necessarily simplify the problem at 
hand, particularly because inversion of a matrix is com- 
putationally intensive and care must be taken in order 
to achieve self-consistency in the eigenvalue. We now 
describe an efficient way to go beyond the standard ap- 
proach to address each of these issues in turn. 

In a plane- wave basis, the inversion of e — can be 
approximated accurately and simply in a way which will 
save considerable computational effort, especially when 
used in conjunction with conjugate gradient techniques. 
Specifically, because the Q-space contains the high mo- 
mentum plane-wave states, the kinetic energy dominates 
HQQ , which we may then approximate to be the kinetic 
energy operator HQQ , whose a,P component is 



H 



QQ _ 
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o,Q/3 — 2'5Q'^a,/3- 



(6) 



Here Sa^p is the Kronecker delta and Q is a reciprocal 
lattice vector, and atomic units have been assumed. 

In principle, the eigenvalues appearing in the left side 
of Equation (|^) must be calculated self-consistently, but 
this leads to difficulties. Each eigenstate sought then sees 
a different Hamiltonian and orthonormality is lost. Al- 
ternately, in practice one employs on the left-hand side a 
constant, approximate e which one hopes to be appropri- 
ate for all desired eigenvalues. For us, neither approach is 
satisfactory. To circumvent these difficulties, we linearize 
C'^ in terms of the eigenvalue 
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)H'''cP+0{-^f. (7) 



Note that in a plane-wave basis such an expansion will 
converge quickly, as the cutoff in the P-space is, in gen- 
eral, much larger than the eigenvalues of interest. 

Substituting Equation (^, to linear order in e, into 
Equation (^, gives an effective generalized eigenvalue 
equation. 



Here 



is the effective single-particle Hamiltonian and 



(8) 



(9) 



(10) 



is a positive-definite overlap matrix. This equation will 
generate errors of order (e/E^)^ in the eigenvalue, where 
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is the cutoff energy for the P-space. We below, show 
that this truncation error proves to be a good measure of 
the error in the energy per atom. Thus, our method pro- 
vides for accurate a priori predictions for transferability 
errors. 

In addition to computing the eigenvalue spectrum, we 
often require access to the eigenstates, particularly for 
density functional theory calculations which require self- 
consistent solution of the electronic states within a po- 
tential dependent upon those states. Calculation of the 
self-consistent potential requires not only the P-space 
components from (^) but also the Q-space compo- 
nents. The most direct choice for generating the C'^ 
is to use Equation (0); however, we wish to eventually 
apply this technique to the conjugate gradient method 
and defining C"' as in Equation (M), no longer preserves 
orthonormality for the full eigenvector set {C} and post- 
reorthonormalization can make the conjugate gradient 
procedure unstable. 

To avoid these difficulties, we note that the solutions 
to the generalize eigenvalue problem (Q) automatically 
satisfy 



a 



where the subscripts indicate the coefficients for individ- 
ual states i and j and we have substituted the definition 
(p^). Regrouping terms, we find that identically 



H. 



so that the set of complete states {C} will be exactly 
orthonormal provided we make the identification, 



(11) 



which is nothing other than (0) truncated at zeroth or- 
der. We thus conclude that to avoid issues of orthonor- 
mality to allow simple implementation of conjugate gra- 
dient techniques, one must construct at this order. 
Accordingly, we employ (^l]) as the construction for 
throughout the reminder of this work. 



B Connections with Lowdin Perturbation Theory 

To carry out full density-functional calculations based 
upon the results from the previous section, one in prin- 
ciple first would solve (H) using direct diagonalization 
techniques to find and then construct C"^ from (11). 
Next, from the wave functions, one would update the 
charge density, Hartree and exchange-correlation poten- 
tials, recompute H^^ , H'^^ ^ and O. Finally, one would 



orytfa ■sdjich was used in the mid 1980's by a number of 
groupstil'LJ with regard to electronic structure calcula- 



iterate this procedure to self-consistency. As this ap- 
proach is quite similar to the Lowdin perturbation the- 

tions, this section briefly reviews the Lowdin approach 
as described in Reference and discuss the major dif- 
ferences between that approach and ours. The following 
section. Section tV, shows how conjugate gradient tech- 
niques may be applied directly to our approach but not 
to Lowdin perturbation theory. 

Lowdin perturbation theory also decomposes behavior 
into high and low momentum plane- wave states. This 
approach also solves an eigenvalue problem 

for a renormalized Hamiltonian U^^ similar to that of 
Equation (S), 



U 
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e-H's 



QQ 



H^P, (12) 



Self-consistency is reached by calculating the charge den- 
sity either with alone, or with both and , with 



^QP^f 



(13) 



Here U'^^ is similar to 11^^, except that the first index 
is in th e Q-s pace. For the reasons discussed above in 
Section III A , if C'^ is employed, then the wave functions 
must be re-orthonormalized. 

In practice,EJ the value employed for e in Equation ( |T^ ) 
depends upon the physics under exploration. Thus, prior 
to the calculation, the shape of the band structure should 
be approximately known and should have a relatively 
narrow band-width, requiring great care for the study 
of systems under unusual physics conditions such as ex- 
tremely high pressures. Moreover, the approximation of 
e by a constant e is always a first order error in the solu- 
tion to the eigenvalue. This can simply be seen by Taylor 
expanding Equations (|l^ ) and ( p^ in a similar fashion 
as Equations (^ and (|^ If e is close to e, then this er- 
ror is small as is the next order. However, if there is a 
large discrepancy between e and e, then these errors are 
manifestly first order. This is particularly troublesome 
for calculations including core electrons where the core 
states energies have large separations from the valence 
states. We thus expect Lowdin perturbation to be useful 
primarily only within a pseudopotential framework. In 
our approach, by contrast, all band energies are treat- 
ing on an equal footing to second order. This is because 
Equation (g) includes all terms through first order by 
shifting some terms to the right-hand side to form the 
generalized eigenvalue problem. 

The great advantage of Lowdin perturbation theory for 
its time is that it reduces the time for standard diagonal- 
ization from {N^ + N'^)^ to (7V^)^. However, in the mid 
1980's standard diagonalization techniques were dropped 
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in lieu of much mnce lefficient Car-ParrinelloiB0 and con- 
jugate gradientEj'ESEj techniques, which, being based on 
minimization, require the eigenvalue equations to repre- 
sent the variational derivative of some energy functional. 
It is not clear that the Hamiltonian used within Lowdin 
perturbation theory is a good approximation to the vari- 
ational derivative of an energy functional. What partic- 
ularly complicates this is the dependence of the Hamilto- 
nian on the eigenvalue and that the unfolding of the wave 
vector into the Q-space does not preserve orthonormality. 
Therefore, it is unclear how one would employ Lowdin 
perturbation theory in conjunction with such optimized 
minimization techniques. The next section demonstrates 
that, in contrast, our approach represents a very good 
approximation to the variational derivative of an energy 
functional and, because it preserves orthonormality, can 
be used directly with optimized minimization techniques. 



IV APPLICATIONS TO CONJUGATE 
GRADIENT TECHNIQUES 

We first discuss applicability to a tra ditional conju- 

and generalize 

to the analytic continued approac" 



and 



gate gradient techniques in SectioiJ^^ 

chESMin Section |IVB. 
Section ^ gives details of our computational implemen- 
tation and provides specific computational studies of the 
various approximations and the convergence of the ap- 
proach. 



A Application to Traditional Conjugate Gradient 
Techniques 

Again, we seek the minimum of an energy functional 
of the form i?[{C}] as a function of a set of wave function 
coefficients {C}. To formulate the separation into low- 
and high- spatial frequency components as a variational 
principle associated with an energy functional, we express 
the energy functional directly in terms of the {C^} and 
{CQ} before taking any variations, E[{C^},{C'^}]. We 
now seek the minimum of this functional subject to the 
orthonormalization condition 



cpocf 



(14) 



Here O is an overlap matrix, which should give the cor- 
rect overlap matrix at the minimum of the energy func- 
tional. Section ^ will describe the overlap matrix that 
we use in our calculations. Finally, the {C^} are left 
unconstrained. 

The sets {C^} and {C^} are treated as independent 
variables. However, at each {C^} point, the functional 
is directly minimized with respect to the set {C^}. The 
gradients are 



dE 



da 



H^PCf + HPQcf - eOCf (15) 



dE 

dcP^ 



(16) 



where the second term is set to zero, emphasizing that 
the energy functional is minimized with respect to {C^} 
at each point {C^}. 

The flowchart in Figure ^ describes the minimiza- 
tion procedure within the traditional conjugate gradient 
framework. Given a point at {C^} and {C*^}, calculate 
the energy. Next, calculate the gradient for the set {C^}, 
Equation (|l|) or (|l|), and hold the set {C"?} fixed. Min- 
imize the energy functional for the set {C^} only, along 
the conjugate direction {X^}. After the functional is 
minimized along this particular direction, the functional 
is then minimized with respect to {C"^}, Equation ([l^), 
by setting 



(17) 



Finally, this process iterates until the energy functional 
reaches its minimum. 



t 

E[{C}] 



SE[(CH 



min{E[{C^+AX'')]) 



t 



FIG. 1: Flow chart for application to traditional conju- 
gate gradient techniques. 



In calculating C"^, we have again approximated H'^'^ 
by H^'^ , and therefore do not exactly minimize the func- 
tional with respect to . Although the form for is 
not that which exactly minimizes our functional, we have 
found the approximation to be sufficiently close that the 
conjugate gradient technique is quite stable and efficient. 
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Using such a form for C"^, the gradient for indeed 
becomes 



dE 
dCT 



= RPPcf + H''ffCf -eOCf, (18) 



where H'^-f^ is as defined in Equation (^. Thus, we find 
that, unUke the Lowdin p ertur bation theory approach, 
the prescription in Section III A represents a set of equa- 
tions that is, to a very good approximation, a variational 
principle and therefore to be amenable to solution with 
conjugate gradient methods. 

To underscore the effectiveness of conjugate-gradient 
methods for use within our framework. Figures ^ and ^ 
show the iterative convergence of the totaLenergy within 
the local density approximation (LDA)l3 to density- 
functional theory for cubic hydrogen and an eight atom 
cell of fee carbon, respectively. The calculations use the 
preconditioner of Teter, Payne and Allen,L3 cutoffs of 
= 8 H and = 16 R for hydrogen and E[' = 30 H 
and E^ = 60 H for carbon, and the overlap operator 
O from Section |^. The hydrogen calculation uses the 
pure|-Goulomb potential and an 8 x 8 x 8 Monkhorst- 
Packe3 sampling grid, and the carbon calculation uses 
the Goedecker, Teter and Hunter pseudopotentiaEj with 
simple F-point sampling. Finally, the cubic lattice con- 
stants were 2.760 bohr and 6.746 bohr for sc hydrogen 
and fee carbon, respectively. 



new approach, with cutoffs stated above, to that with the 
traditional approach, with cutoff E*^'^'^ = E^ . The fully 
converged results are considered to be that of the tra- 
ditional approach with energy cutoff E^°™ = E^ . The 
results for the new procedure show a vast improvement. 
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FIG. 3: Energy error versus iteration number for the new 
method, within the conjugate gradient framework for fee 
carbon when using a nonlocal pseudopotential. 

As noted above, one of the advantages of our approach 
is that it allows for an a priori estimate of the transfer- 
ability errors. We define this transferability prediction as 
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FIG. 2: Energy error versus iteration number for the 
new method, within the conjugate gradient framework 
for simple cubic hydrogen when using the Coulomb po- 
tential. The first few data points converge slowly in order 
to accurately converge the fillings. 

It can be seen that the conjugate gradient technique is 
stable and that fluctuations occur only in the final stages 
when the error is AE/Emin ~ 10"^, well within accept- 
able convergence. These fluctuations are associated with 
the fact that H^'^ is used instead of H^^. In general, 
we have found that the stability of the conjugate gradi- 
ent does not depend on the approximation for C"^, but 
mainly depends on how and when we update . Ta- 
ble I compares the error in the energy per atom with the 



^ ^ ^states 
states 



{estates/E^y, (19) 



where the sum is over occupied states. For the spe- 
cific cases of the present two calculations, our prediction 
gives AEP^^J^"* « 0.5(0.5/8)2 ^ 2mH for hydrogen and 
AEfH^^"* ^2x 0.5(0.5/30)2 ^ o.3TOi? for carbon, where 
for simplicity we consider only the 2 s-states. These pre- 
dictions also appear in Table |[ The predicted values are 
very similar to AEnew, demonstrating that this predic- 
tion, which may be applied in any environment, gives 
sensible estimates of the errors. 





AEtradiEc = E':) 




A rppredtcted 


Hydrogen (sc) 


9niH 


ImH 


2mH 


Carbon (fee) 


l&mH 


0.3mH 


0.3m// 



TABLE I: The error in calculating the energy per atom 
for simple cubic hydrogen and fee carbon. The first col- 
umn shows the error in energy when using the traditional 
approach. The next column shows the error in energy 
when using the new approach. The final column shows 
the transferability prediction of the new approach. Equa- 
tion (|l|). 



Sections V C and VI show that in practice the present 
method is particularly beneficial for aH-electron calcula- 
tions of first row elements, where the reduction in the 
cutoffs is much more dramatic than the cases presented 
here. 
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B Application in the Analytic Continued Approach 



Finally, the solution for becomes 



Our procedure can be easily incorporated into the an- 
alytic continued conjugate gradient approach^ In such 
an approach, constrained variables, {C}, are not mini- 
mized, but unconstrained variables, {Y}, are minimized. 
This approach allows for much better convergence, as all 
directions are allowed in the search space when minimiz- 
ing the energy functional.E^ 

In the analytic continued approach, the sets {C} and 
{Y} are related by 



C = Yu-iv\ 



(20) 



where bold faced quantities refer to the expansion co- 
efficients for all states gathered into matrices, each of 
whose columns represents a particular state. Here, V is 
a unitaHV-feransformation which allows for subspace ro- 
tationsE3'Ell and u is the expectation value of Y with 
respect to the overlap matrix O, 



= Y^OY. 



(21) 



Note that Equation (gO|) is constructed so that the C 
automatically satisfy the orthogonality constraint, with 
regard to O, as direct substitution verifies. This ap- 
proach allows the total energy to be found by minimizing 
the energy E[Y] directly using standard, unconstrained 
conjugate-gradient techniques. 
In our case, the gradient becomes 



dE 

gypi 



(1 - OC^C^^)(//^^ + H'ff)C^FVU^ 
+ OCPVQ{V'^[H,F]V), (22) 



where H'^^^ is as defined in Equation (^), F is a diagonal 
"filling" matrix composed of the state occupancies, O is 
the overlap matrix, defined-|Si«iilarly to Equation (10), 
and QO is the (5-operator,E£lEZl defined by 



{W^Q{A)W)^ 



{W'^AW)^ 



where W is the unitary, column order, eigenvector ma- 
trix for u and /i are the corresponding eigenvalues. In 
Equation (|2^), the Q-operator operates on the matrix 
V'^[H, F]V , where [i?, F] is the commutator between the 
filling matrix F and the subspace Hamiltonian matrix. 



Now, u is redefined as 



H 



•ff\ 



yP^OY' 



(23) 



and is related to by 



C^=Y^u-5V^t. 



(24) 



We note that H"^^^ in Equation ( p2[ ) acts as a pseu- 
dopotential generated directly from the crystal environ- 
ment and not from a referenced atomic calculation which 
one must hope to be transferable. The overlap matrix O 
acts to preserve the correct orthonormalization for the 
P-space wave function. This overlap matrix is similar 
to those found in the Ultra-Soft Pseudopotentialsc^ and 
Projected Augmented WaveE3 methods. In those meth- 
ods, the softening procedure is generated in a real space 
formalism, from a specific reference state, and depends on 
the choice of a cutoff in real space cutoff. In this work, 
the softening is done in a momentum space formalism 
without artificial core radii, making the procedure ideal 
for high-pressure studies where the loss of distinction be- 
tween core and valence regions is precisely the physics of 
interest. As a final advantage of the present approach. 
Equation ( p^ ) gives a direct measure of transferability 
errors. 



V COMPUTATIONAL DETAILS 



We will now describe briefly efficient calculation of spe- 
cific terms in our procedure. The text focuses on all- 
electron calculations as these are of interest in the present 
work. The appendix discusses details for the use of this 
approach with hard but transferable Kleinman-Bylander 
pseudopotentials . 



A Form for Overlap Matrix 



The overlap matrix has the form 



O = 1 



H 



PQ 



H. 



pq\ 



1 



1 



QP 

h,xc 



H. 



QP\ 



(25) 

where -j^^ is the inverse of the kinetic energy, Hion 
is the contribution from the local ionic potential and 
Hh.xc is the contribution from the Hartree and exchange- 
correlation potentials, which depends on the density. 
There are a number of viable options for implementing 
O. We have explored both the option of fixing Hh,xc: 
and thus O to some value and the option of updating it 
self-consistently. 

The benefit of holding the overlap matrix fixed is that 



do 



= 0, 



(26) 



so that Equation ([1^) is the exact variational derivative 
of the energy functional. The appendix further shows 
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that such a fixed operator makes efficient implementa- 
tion of Kleinman-Bylander pseudopotentials possible. A 
clear choice for a fixed overlap matrix is to hold Hh.xc in 
Equation ( ^5|) to that of the "free-atom" crystal, a crys- 
tal whose charge density is just the superposition of the 
charge densities from isolated atoms. In such a case, the 
overlap matrix becomes 



where H'^ 



1 



1 



jjQQ jjQQ 



jjcr-aUQP^ (27) 



-at,PQ _ TTCr-at,PQ 



IS 



Hartree and exchange-correlation potential for the "free- 
atom" crystal. One possible cause for concern when us- 
ing a fixed overlap matrix is that the orthonormalization 
is no longer exact (so that, for instance, the electronic 
charge will not be precisely conserved) because H'^^ in 
Equation (24) will differ from its free-atom crystal value. 
However, we have found this to never affect the conver- 
gence of the calculations and to yield quite accurate re- 
sults, as Figures || and |^ and Table || show. 

Alternately, exact orthonormalization can be preserved 
when using a fixed overlap operator if one chooses to re- 
define as 



1 



H, 



QQ 



jjcr-at.QPQP 



(28) 



While this preserves orthonormality among the states, it 
does give a slightly different approximation to C^, which 
causes a slight change in the H'^^-^ which should be used 
in (^. We find, however, that ignoring this discrepancy 
makes no practical difference in the final results. With 
this alternate approach, the convergence is quite similar 
to those of Figures |^ and and the energy per atom is 
within 1/i H, as compared to using Equation (M for C<3. 
Another benefit of using the form Equation (|28|) comes 
when using non-local pseudopotentials, as this form is 
computationally more efficient. (See the appendix for 
details.) 

Finally, the H'^f^ which should be used with the con- 
struction in Equation (E 



We have found somewhat to our surprise, however, 
that updating the overlap matrix self-consistently while 
completely ignoring these contributions to the gradient, 
that the gradient technique remains stable and maintains 
its good convergence, as Figures ^ and || show. In these 
calculations, H'j^~^* was used only at the first iteration to 
initialize the calculation. After that, the overlap was up- 
dated self-consistently and the gradient was always cal- 
culated ignoring the aforementioned contributions from 
the changes in O. The convergence using this technique 
is quite similar to that which holds O fixed. The only no- 
ticeable difference occurs in the asymptotic region, where 
now slightly larger ffuctuations exist as ffuctuations in 
the overlap matrix are not properly taken into account. 
Finally, we find the final energy per atom again to be 
within 1/xH of all of the approaches above. 
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FIG. 4: Energy convergence versus iteration for simple 
cubic hydrogen when the overlap matrix is updated self- 
consistently, but its variational derivative is not taken 
into account. 



IS 



H-ff = _hPQ 



-at.QP 



(29) 



Maintaining this consistency comes at the cost of intro- 
ducing a slightly non-Hermitian Hamiltonian. Despite 
this, we again have not experienced any practical prob- 
lems working with this form. 

The second option for implementing the overlap ma- 
trix is to use the Hartree and exchange correlation po- 
tentials self-consistently. When doing so, the gradient of 
the total energy involves terms including the derivative 
of the overlap matrix. Unfortunately, we have yet to find 
a method for calculating these contributions to the gra- 
dient which does not require multiplication by matrices 
of size X , a completely impractical feat for all 
but the most trivial of systems. 
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FIG. 5: Energy convergence versus iteration for fee car- 
bon when the overlap matrix is updated self-consistently, 
but its variational derivative is not taken into account. 
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B Numerical Implementation 

Specifying the mathematical form of the overlap ma- 
trix in one of the forms from the previous section and 
thus also the forms for and H*^^^ leaves the task of 
evaluating these and related quantities numerically. This 
section first discusses efficient methods for this evaluation 
and then describes how to build conjugate gradient min- 
imization from the ability to perform these evaluations. 

The term OY^ must be calculated for defining m, 
Equation (p^), and in the gradient, Equation (p^). The 
most efficient method for calculating this term is to per- 
form four Fourier transforms per electronic wave func- 
tion, always multiply terms in the space in which they 
are diagonal. We find it useful to keep this term in stor- 
age so as to only calculate it once in a conjugate gradient 
loop. Such storage is minimal, as will be seen. 

To reduce the time spent in Fourier transforms, we 
make the approximation of using a slightly smaller 
Fourier transform grid than might normally be used. 
The gradient, the high momentum coefficients C^, and 
the overlap matrix only involve momentum transfers of 
P + Q. Higher momentum transfers do not occur be- 
cause is approximated by H^'^, which is diagonal. 
Accordingly, rather than Fourier grids of size 2Q, we em- 
ploy grids of size P + Q. We find that using the full 
Fourier transform grid decreases the energy minimally: 
for hydrogen and carbon doing so decreases the energy 
per atom by only 0.5/.tH and O.SfiH, respectively. 

The reduction of the Fourier transform grid to size 
P + Q from a size of 2Q, which would be needed in a 
non-perturbative approach, saves significant storage in 
the ratio 



FFT„ 



FFT- 



full 




Ep 
Eq 



El 

Eq 



(30) 

In our boron calculation, this is a factor of nearly three, 
which can make the difference between being able to per- 
form the calculation on a workstation or requiring paral- 
lel computing capability. 

We find the following procedure to be very stable for 
performing conjugate gradient minimization. First, given 
generate according to Equation ( p4[ ) or (28). 
Next, calculate the total energy and the gradient for Y-'^ 
at this point from Equation (^). Holding fixed, min- 
imize along the search direction, and generate a new 
Y''^, and iterate to convergence. 

This procedure for solving for the P-space and gen- 
erating the Q-space, when needed, is very beneficial for 
plane-wave calculations when using either the Coulomb 
potential or a very hard, yet highly transferable, pseu- 
dopotentials because, generally for such calculations the 
computational bottleneck is memory and this approach 
requires permanent storage only in the P-space. (The 
Q-space wave functions for a particular band can be gen- 
erated only when needed.) Moreover, in terms of com- 



putational time, a savings of 1 -t- Nq/Np for all direct 
matrix multiplications occurs, and the P-space approach 
is therefore also more efficient in terms of computational 
time for systems with many bands. 



C Application to the Coulomb Potential 

To demonstrate the benefits of our approach, this sec- 
tion presents all-electron calculations for fee boron at two 
different volumes {Vi = 3.83 and V2 = 3.5 A^), ex- 
ploring convergence with energy cutoff, usage of memory 
and computational time, and our ability to predict trans- 
ferability errors. We carry out all calculations within the 
local density approximation with ten special k-points in 
the reduced Brillouin zone and a Fermi surface integra- 
tion temperature of kT « 0.0037 H. 
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FIG. 6: Energy per atom as a function of cutoff energy for 
fee boron, when using the Coulomb potential in the tra- 
ditional plane-wave approach. The squares correspond 
Vi = 3.83 A^ and the circles correspond to V2 = 3.5 A^. 

Figure |^ shows the convergence of the energy as a 
function of cutoff for the direct plane wave approach. 
The data in the figure demonstrate that the traditional 
plane- wave approach converges quite slowly, requiring an 
energy cutoff of over 12000 Ryd to bring the total energy 
to within 0.1 eV per atom. However, as is well known, 
the higher-energy plane wave states provide convergence 
mostly to the inert region of the core near the nucleus 
and thus physically meaningful energy differences con- 
verge much more rapidly than the absolute total energy. 
As a more meaningful reference, Figure ^ shows the con- 
vergence of the energy difference between volume states 
Vi and V2 using the direct plane wave approach. 

Figure ^ shows the convergence with P-space cutoff of 
the total energy when using the perturbative approach 
with a Q-space cutoff of 7000 Ryd and employing H'j^^~°'* 
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for construction of O, and H'^f^ (Equation (p9|)). 
The figure shows that there is almost no difference in the 
total energy predictions (circles and squares for the two 
volume states) when reducing the cutoff in the P-space 
to 3500 Ryd and that reducing the P-space cutoff down 
to 1000 Ryd gives just as accurate total energies as a 
6000 Ryd cutoff in the traditional approach. The figure 
also shows that the transferability prediction based on 
the two core electrons (diamonds) gives the correct order 
of magnitude and is never off by more than a factor of 
three in this case. 
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FIG. 7: Energy difference as a function of cutoff for fee 
boron, when using the traditional plane-wave approach. 
The volumes per atom are Vi = 3.83 and V2 = 3.5 A"^. 



2.7 which largely compensates the need for a few more 
Fourier transforms in generating and OY^. 
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FIG. 8: Energy convergence versus P-space cutoff, using 
the new method, when = 7000 Ryd. The squares 
correspond to Vi = 3.83 A"^, the circles correspond to 
V2 = 3.5 A'^ and the diamonds correspond to the trans- 



ferability prediction. Equation ( |19| ) . The horizontal lines 
correspond to the energy calculated using the traditional 
approach at a cutoff energy of 7000 Ryd. 



Figure ^ shows the performance of the new approach in 
computing energy differences. Energy differences in the 
perturbative approach tend to fluctuate slightly and lack 
the monotonic behavior which the traditional approach 
exhibits due to loss of ability to describe physical changes 
in the high momentum states. Perturbation theory re- 
covers these changes nearly perfectly and so the resulting 
error in energy difference is smaller and loses systematic 
behavior. From these results we conclude that a cutoff in 
the P-space of 1200 Ryd and of 7000 Ryd in the Q-space 
is more than sufficient to predict accurate energy differ- 
ences and to capture any unforeseen physical processes 
which come into play. 

To make a direct comparison of computational sav- 
ings, we note that the total energy and energy difference, 
when using the new method with = 1200 Ryd and 
E^ = 7000 Ryd, is very close to the traditional approach 
with energy cutoff El^""^ = 6500 Ryd. At these cutoffs, 
the memory savings from the perturbative approach is 
a factor of 12.5 for the wave functions and 2.7 for the 
FFT's. The time for matrix multiplications decreases by 
a factor of 5.4 and for Fourier transforms by a factor of 
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FIG. 9: Energy difference in the new approach as a func- 
tion of the P-space cutoff. 
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VI APPLICATION TO HIGH PRESSURE 
BORON 



As a physical application, we study the high-pressure 
phases of boron, whose rich physics remains mysterious. 
Boron is known to exhibit a semiconducting to metallic 
phase transition at w 170 GPaJj and under high pressure 
and low temperatures, boron superconducts.BI Previous 
theoretical studies fit boron under pressure include that 
of Mailhiot et. al.^ who used both the LMTO method 
and the pseudopotential method. They have found their 
pseudopotential results to be more reliable and predict a 
sequence of structural phase transitions with increasing 
pressure from the icosahedral structure (ai2-B) to a body 
centered tetragonal structure (bet) to the face-centered 
cubic structure (fee). They find the ai2-B bet tran- 
sition to occur at 210 GPa, and established this as an 
upper bound for the semiconducting to metallic phase 
transition, within the local density approximation. 

To calculate the energy for different phases in boron 
we use our perturbative potential approach with the an- 
alytically continued conjugate gradient approach as de- 
scribed in Section IV with O, and H'^^^ calculated 



as in Section V C , All calculations employ the local den- 
sity appmximation with the parameterization of Perdew- 
ZungerESl to the Ceperley-Aldcrcd exchange-correlation 
energy. 

WCi^lst) consider other structures than explored previ- 
ouslyJj'tHl not only fee, bet, and q;i2-B structures but also 
orthorhombic structures of boron with the a-gallium and 
/3-gallium basis, which we denote as a^a-B and /3ga-B, 
respectively. We consider these orthorhorohic structures 
for a number of reasons. Zhu and HenleyCj have shown 
that for 6-coordinated boron under high pressure the 
triangular lattice is a good candidate for a low energy, 
high pressure phase, within the-density functional the- 
ory pseudopotential ixamejaaprkcJ and furthermore such 
sheets tend to buckle 0'c£Il3 Because the gallium struc- 
ture can be viewed as stacked sheets of a highly buckled 
triangular lattice, such a structure is a good candidate 
for a high pressure phase of boron. Secondly, gallium is 
another Group III element which has as its ground-state 
structure a-gaUium at low pressures and /3-gallium at 
high pressures ;Eij and being under boron in the periodic 
table supports the notion that such structures are good 
candidates for a high pressure phase of boron. Below, we 
show that these structures play quite important roles. 

For the bet structure, we use the same c/a — 0.65 ratio 
as was used in Reference ||, where this ratio was shown 
to give the ground state at a volume « 3.91 within 
the pseudopotential calculation. We have tested this ra- 
tio for a number of high pressure states and have found 
the ideal ratio to remain constant over a wide range of 
pressures. Hence, we hold c/a at this value for all cal- 
culations in our study, and thus our bet energies results 
represent a quite close upper-bound. 

For the q;i2-B structure, we use the structure listed 



in Table Suof Reference |T^, the same used by Mail- 
hiot et. ala This structure is rhombohedral with angle 
a ~ 58.06°. The basis locations for the atomic positions 
are 

±{xixizi; xizixi, zixixi] X2X2Z2] X2Z2X2] Z2X2X2), 

where xi = 0.0104, zi = -0.3427, X2 = 0.2206 and Z2 = 
—0.3677. The space group for the q;i2-B structure is 
R3m. Note that we may regard this q;i2-B structure as 
representative of the low-pressure states of boron which 
are all governed by the icosahedron. 

For the orthorhombic structures, we^se the Gallium 
structures as found in Wyckoff 's book.c3 For the aga-B 
structure, we hold the lattice vector ratios to a/b = 
0.99867 and a/c = 0.590035 and place the atomic ba- 
sis at coordinates 

±{0,u,v; l/2,u+l/2,w; 1/2, u, 1/2; 0, u+1/2, 1/2-w), 

where u — 0.0785 and v — 0.1525. The space group is 
Cmca. For the /3ga-B structure, we hold the the lattice 
vector ratios at a/b = 0.3567 and a/c= 0.9148 and place 
the atomic basis at 

±(0,u,l/4; 1/2,^^+1/2,1/4), 

where u = 0.133. The space group is Cmcm. 

Note that we do not optimize the lattice vector ra- 
tios or the internal coordinates for the orthorhombic 
structures for boron but here hold them fixed to that 
of gallium at standard conditions. Therefore, the ener- 
gies which we report for these structures should be re- 
garded simply as upper bounds. Similarly, we do not 
optimize the lattice vectors and/or internal coordinates 
for the bet and ai2-B structures, however this—latter 
optimization proves to be minor. Zhao and LuEHI have 
recently performed pseudopotential calculations, for the 
fee, bet and ai2-B structures, and found that fully re- 
laxing these structures only decrease the transition pres- 
sure (bet ai2-B) by 5%, when using the local density 
approximation, as compared to Reference ^. Therefore, 
because these latter effects are minor near the semicon- 
ducting/metallic transition pressure and in order to allow 
our method a greater comparison with pseudopotential 
results we choose to hold the parameters for the bet and 
Q!i2-B structures fixed to those used in Reference ^. 



A Convergence 

We now establish plane wave cutoffs for our specific 
study which will both make for practical calculations 
(particularly for the ai2-B structure) and give accurate 
structural energy differences. For this we shall compare 
the energy per atom as a function of volume for the fee, 
bet and ai2-B. As these are convergence tests, we employ 
relatively modest Brillouin sampling: 4x4x4, 4x4x6, 
and 1x1x1 (the F point) meshes for the fee, bet and 
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ai2-B, respectively. To facilitate integration over the 
Fermi surface, we employ an electronic temperature of 
kT w 0.0037 H. 

Figure summarizes the results for the three struc- 



tures when using the plane wave cutoffs from Section V C , 
EQ = 7000 Ryd and = 1200 Ryd. The curves from 
the figure are simple fits to 



> 

s 

o 
a 

a 



(31) 





• fee (1200/7000) 

• bet (1200/7000) 

4 a1 2-6(1200/7000) 


\ \ 




■,\ if 




% 








^ 



> 

CD 



o 

< 

CD 

a 

CD 



"A Y 


« fee (200/1200) 
. bet (200/1 200) 
« a1 2-6 (200/1 200) , 

--■ fee (1200/7000) 
bot (1200/7000) 

a1 2-6(1 200/7000), 


-A • 

'. \ Y 

A 





4 5 6 7 

Volume per Atom [A^] 



4 5 6 7 

Volume per Atom [A^] 



FIG. 11: Energy per atom versus volume per atom. 
The results from new method using the cutoffs £;f =200 
Ryd and i?f^=12GG Ryd are displayed as circles (fee), 
stars (bet) and diamonds (ai2-B). The results from new 
method using the cutoffs E^f =1200 Ryd and £;^=700G 
Ryd are displayed as lines by fitting the data to Equa- 
tion (|3ll): fee (dashed), bet (dotted), ai2-B (solid). 



FIG. 10: Energy per atom versus volume per atom for 
the fee (circles), bet (stars) and ai2-B (diamonds) struc- 
tures. Energy cutoff for the P-space is 1200 Ryd and for 
the Q-space 7000 Ryd. The curves are fitted to Equa- 
tion 

Having established previously that the above cutoffs 
give extremely good energy differences, we then searched 
for possible reductions in the cutoffs which continue to 
reproduce these results accurately. We find reducing the 
cutoffs to 200 Ryd for the F-space and 1200 Ryd for the 
Q-space to both lead to efficient calculations and to give 
highly accurate results. Figure ^ compares the results 
at this reduced cutoff (data points) with the previous 
results at the higher cutoffs (curves). 

As a comparison, we also calculate the energy when 
using a traditional, direct plane-wave approach with a 
cutoff of 200 Ryd. Figure |l2| shows a similar plot, where 
the points are the plane wave data and the curves rep- 
resent the fit to the fully converged results. Note that, 
although fee structure is fairly well converged at this low 
cutoff, there are rather large errors for the ai2-B struc- 
ture. Thus, such a low cutoff is unreliable for the study 
of such systems. The fact the low cutoff represents one 
structure well and not another demonstrates the unpre- 
dictability of transferability errors, underscores the need 
to predict those errors reliably, and places pseudopoten- 
tial studies of such systems, which assume that errors in 
processes at such high energy scales cancel, in doubt. 
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FIG. 12: Energy per atom versus volume per atom. The 
results from using the traditional approach with cutoff 
E=200 Ryd are displayed as circles (fee), stars (bet) and 
diamonds (ai2-B). The results from new method using 
the cutoffs £:f =1200 Ryd and £:^=7000 Ryd are dis- 
played as lines by fitting the data to Equation (31): fee 
(dashed), bet (dotted), ai2-B (sohd). 



Finally, we note that the memory savings of using the 
new method, with cutoffs of 200 Ryd and 1200 Ryd, ver- 
sus the traditional approach, which would require a cut- 
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off of 1200 Ryd to produce results of similar reliability, 
is quite substantial, a factor of « 15 for wave functions 
and « 3 for the FFT box. The time for the matrix mul- 
tiplication correspondingly reduces by a factor of 6 and 
for the Four ier transforms by a factor of 3, similar to 
Section V C , Because of these savings, particularly in 
memory, all of our calculations below were possible on a 
single desktop computer and there was no need for par- 
allel supercomputing. 



B Results 

Having found appropriate plane wave cutoffs, we next 
converged both the Brillouin sampling and the fictitious 
Fermi temperature. We find that a fictitious Fermi tem- 
perature of kT = O.OOli? converges the total energy to a 
within a few tenths of a millihartree per atom {Q.2mH for 
the fee structure). For the fee and bet structures, Bril- 
louin sampling on 6x6x6 and 4x4x6 meshes, respec- 
tively, suffice to converge the total energy per atom to 
within the same tolerance. We then generate meshes for 
the remaining structures by maintaining the same recip- 
rocal space sampling density as closely as possible, re- 
sulting in meshes of size 2x2x2, 4x4x4 and 6x4x6 for 
the ai2-B, a^a-B and /3ga-B structures, respectively. 



> 



< 




o fee 

. bet 

<f a12-B 

< aga-B 

► Pga-B 



3 4 5 6 7 8 

Volume per Atom [A^] 

FIG. 13: Energy per atom versus volume per atom, for 
the fee (circles), bet (stars), ai2-B (diamond), a^a-B (left 
facing triangles) and /?ga-B (right facing triangles) struc- 
tures. The results are from the new method using the 
cutoffs EP=2QQ Ryd and E^^Um Ryd. 

Figure |l^ shows our final results, which are converged 
to within a few tenth of a millihartree per atom. Most no- 
tably, we find a new phase, the a^a-B , to be the lowest in 
energy phase under high pressure. a^a-B boron becomes 
energetically favorable over ai2-B boron at a pressure 



of 71 GPa, and remains the energetically favored phase 
up to the highest pressures we have studied, which are 
over 500 GPa, where the pressure is determined from the 
enthalpy. Although there are other stable structures gov- 
erned by the icosahedron in boron, due to the dramatic 
rise in energy of the ai2-B structure, we do not expect 
these alternate structures to hpj-e significantly lower en- 
ergies. Moreover, Zhao and Lull3 have preformed prelim- 
inary pseudopotential calculations with the /3-Bio5 unit 
cell, part of the icosahedral family, and find similar trends 
between this structure and the ai2-B structure. Finally, 
optimizing the lattice vectors ratios for the orthorhombic 
phases will only make them more favorable. We, there- 
fore, expect our prediction of a phase transition from the 
icosahedral family to the a^a-B structure to be robust. 

As a iwint of comparison to the pseudopotential cal- 
culation,Ll which had not considered the orthorhombic 
structures, ignoring those structures, we also find the se- 
quence of transitions ai2-B bet — > fee. Our value of 
194 GPa for the first of these transitions agrees reason- 
ably well with the pseudopotential result of 210 GPa. 
However, our prediction of « 475 GPa for the second 
transition differs substantial from the pseudopotential re- 
sult of 360 GPa. These discrepancy most likely results 
from the inadequacies of the pseudopotential to effec- 
tively account for the response within the core. The re- 
sults in Figure |l^ demonstrates the importance of mo- 
mentum transfers over 200 Ryd and hence variations < 
0.44 ao, where ao is the Bohr radius, well within the core 
region. This underscores the importance of control over 
transferability and properly accounting for the core re- 
gion when studying high pressure systems. 
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FIG. 14: Density of states for the valence electrons for 
the ttga-B structure. The straight line corresponds to the 
chemical potential at kT — G.OOli?. 
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FIG. 15: Density of states for the valence electrons for 
the /?gQ-B structure. The straight line corresponds to the 
chemical potential at kT = O.OOli?. 

Given tlie availability of high-pressure resistivity data 
for boronp it is interesting to consider the electronic 
structure of the hitherto unexplored ctga-B phases for 
boron. Figures |l^ and |l5| show low-resolution plots of the 
density of states for a^a-B and /3ga-B boron at 250 GPa 
and 217 GPa, respectively, as generated using the Bril- 
louin zone sampling from the total energy calculations. 
The plots have low resolution because, as is well known, 
meshes which give reliable total energies often are too 
modest to give detailed density of states. (Providing 
smooth curves required a Gaussian broadening of width 
1 eV.) The figures clearly show that, whereas the /3ga-B 
phase is semimetallic, the density of states of the a^a-B 
phase clearly exhibits a gap at the Fermi level, even at 
high pressures. 

Experimentally, the room-temperature resistance of 
boron decreases discontinuously as a fmiction of pres- 
sure at 30 GPa, 110 GPa and 170 GPa.B In comparing 
our results to these data, it is important to recall that 
the ai2-B structure is only meant to be representative 
of the icosahedral family of low-pressure structures and 
that, with optimization of the aspect ratios and internal 
coordinates, the curves for the a^a-B and /3ga-B struc- 
tures will displace to even lower energies. In particular, 
a downward shift of only 0.15 eV ( 10%) in the aga-B 
curve would lower our prediction for the ai2-B-^-ag£i-B 
transition from 71 GPa to 30 GPa. Thus, while further 
investigation is clearly needed, these data suggest that a 
structural transition from the icosahedral family to the 
orthorhombic phase is a viable candidate for the resis- 
tance discontinuity observed in boron at 30 GPa. It is 
also interesting to note that the aga-B structure (Cmca 
space group) was also observed to be stable low energy, 
high pressure phases for both lithium and sodium Jj'cl 
within density functional theory. 



Because boron is observed experimentally to be metal- 
lic above 160 GPa, our results imply that another, yet 
unexplored, structure must become energetically favored 
at this pressure and be responsible for the highest of the 
observed discontinuities in the room-temperature resis- 
tance. Whatever this metallic structure is, it cannot be 
simple as the fee, bet, bcoJ and ideal hcpO structures 
because none are competitive at high pressures. It is 
suggestive, however, that our energetic upper-bound for 
the metallic /3ga-B structure is quite competitive with 
the fee and bet structures at the very highest pressures. 
Although further calculations would certainly be needed 
to verify this conjecture, it is quite possible that the ob- 
served discontinuity in the resistivity at 170 GPa rep- 
resents a structural transition from the semiconducting 
Q!ga-B to the metallic /3ga-B phase. Moreover, it is un- 
clear if pseudopotential calculations will give the accu- 
racy needed in this transition region. 



VII CONCLUSIONS 



We have developed a new method which allows for dra- 
matic reduction in the number of plane waves needed 
in all-electron density functional theory calculations or 
in calculations with hard, but highly transferable pseu- 
dopotentials. The foundation of the approach is to 
treat the higher energy plane wave components of the 
electronic states implicitly through a perturbative ap- 
proach, thereby allowing these components to be recov- 
ered quickly only when needed, thus alleviating mem- 
ory storage requirements by factors as high as fifteen in 
practical calculations. The approach lends itself well to 
current optimized minimization techniques such as con- 
jugate gradient methods and has the benefit of allowing 
quantitative estimates of transferability to new systems 
for which there is little or no experience and for which the 
transferability of pseudopotentials is in question, such as 
condensed matter at high pressure where effects in the 
core region become much more important. 

In addition, the first application of this approach has 
led to new intriguing conjectures as to the structure of 
boron at high pressures. We have carried out the first 
ab initio calculations of the orthorhombic structures of 
boron, which we show to play an important role in the 
high pressure behavior of this material. We find that, 
prior to becoming metallic, boron makes a phase transi- 
tion from the icosahedral family to the semiconducting 
ckga-B structure. The low energy of this semiconducting 
phase indicates that the structure of the superconduct- 
ing phase for boron is more complicated that the simple 
monotonic lattices explored to date. While further cal- 
culations including optimization of the lattice ratios for 
the orthorhombic phases are needed, our results lead to 
the conjecture that the metallic /3ga-B structure is a good 
candidate for the superconducting phase. 
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APPLICATION TO KLEINMAN-BYLANDER 
PSEUDOPOTENTIALS 



This appendix reviews briefly the application of our 
approach to norm-conserving pseudopotentials of the 
Kleinman-Bylander form, where the approach essentially 
becomes a softening procedure without transferability bi- 
ases. The softening comes at the cost of generating a 
generalized eigenvalue problem, siejiilar in spirit to the 
ultrasoft pseudopotefttial (USSP)Eil and projector aug- 
mented wave (PAW)E3 methods. However, unlike the 
USSP and PAW approaches, where the softening takes 
place referenced to a spherical atomic environment, in 
our approach the softening occurs in the full crystalline 
environment and thereby includes nonlocal interactions 
from multiple scattering events from different atoms in 
the crystal but at minimal extra computational cost. 

Here, we develop the approach to work with standard 
separable non-local potentials of the form 

Hr.l{G,G') = \Vat,priG))fpr{Vat,pr{G% (1) 

at^pr 

where the sum is over all atoms in the cell and all projec- 
tors on a given atom. For such potentials it proves much 
more computationally efficient to use the forms similar 



to Equations (t 
O, and H^^ 
takes the form 



, (|2^) and il 
respectively. 



3|) for the overlap matrix 
The overlap matrix then 



O 



1 



the remaining terms contract contract onto the P-space 
only. We now consider each of these in turn. 

The p, jf component of the nonlocal-nonlocal term is 

Onl,nliP,p') ^ ^ \Vat,prip))fprA^\l,p, fpr'{Vat',pr'{p')\, 

(3) 



where 



2^ 



|Kt',pr'(g)). (4) 



This term quantifies how the high-momentum states 
components couple to the ions, including not only cou- 
plings to individual atoms but also to pairs of atoms 
through A p, , thereby including environmental effects 

directly into the potential and enhancing transferability. 
In terms of implementation, this term requires almost 
no additional memory as the only additional storage is 
for A^^^ and scales as iVf, x iVf,, where is the number 
of bands. The computational demands are quite light 
as requires an amount of computing comparable to 
a single orthonormalization of the electronic bands and 
only need be compute once at the very beginning of the 
calculation. 

The local-nonlocal and nonlocal-local terms are just 
hermitian conjugates of each other. Thep, p' component 
of the nonlocal-local term has the form 



where 



\q^^). (6) 



can also be simply calculated in the begin- 
ning of the calculation through Fourier transforms of 

(W(90|(tV)% 

, . The effective Hamiltonian H'^^^ can also be contracted 
^ ■onto the P-space. H'^^^ has the form 



PQ 



cr-at,PQ 
loc 



jjQQjjQQ 



H 



QP 

nl 



H. 



-at,QP 



loc 



H-ff 



H 



PQ 

loc 



1 



where H^^^"" is the local term for the "free-atom" crys- 
tal, a crystal whose charge density is a direct superpo- 
sition of the charge densities from isolated atoms. This 
term includes the local part of the ionic potential, the 
Hartree term and the exchange-correlation term. Ignor- 
ing the identity, expansion of Equation (H) results in four 
terms, which we refer to as nonlocal-nonlocal, nonlocal- 
local, local-nonlocal and local-local, respectively. The 
local-local term is exactly the same as that of the 
Coulomb potential and therefore should be treated as 
in the manuscript. The advantage of employing a fixed, 
rather than self-consistent, overlap matrix is that each of 



loc 



-at,QP 



h: 



nl 



-at,QP 



The nonlocal term on the left-hand side can be contracted 
onto f{':r-at,QP similar fashion as was done to O, 
giving two terms: 



Hnl,nlip,p)^ Y \Vat,pr{p))fprA^llp,fj 



at ' , pr ' 



where 



9 



'{Vat',pr'{p)\, 

(7) 

{Vat,prm ( ) \Vat',pr'{q)), (8) 

2^ 
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and 



at^pr 



where 



{Fl^lrm = {Vat,A€)\ f T^) Hzr\<i- f)- (10) 

V 2^ / 

This will reduce almost all matrix multiplications onto 
the P-space. The only remaining Q-space multiplication 
would be in -ffjoc.n; A final simplification that can be 
made is to take 

rre// TTcr—at^PQ ^ Trcr—at^QP 

and thereby make the full eflfective Hamiltonian con- 
tractible onto the P-space. 
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